#include "BVH.h"
#include "Ray.h"
#include "Console.h"
#include "Scene.h"

#ifdef SSE
#include "SSEStuff.h"
#endif


BVH::BVH(int depth)
: m_depth(depth), m_objects(0)
{
}

BVH::~BVH()
{
	clean();
}

void BVH::clean()
{
	// clean up object list
	int n = m_objects->size();
	for (int i=0; i<n; ++i)
	{
		Object *obj = (*m_objects)[i];
		// if obj is a BVH
		if (obj->getType() == TYPE_BVH)
			delete obj;
	}
	delete m_objects;
	m_objects = 0;
}


// construct the bounding volume hierarchy
void BVH::build(Objects * objs)
{
	// if BVH has already been built, unbuild it
	if (m_objects)
		clean();
	m_objects = new Objects;

	// bounds of BVH are tightly fitting box around all child objects
	getBounds(objs, &m_min, &m_max);

	int n = objs->size();
	if (m_depth < BVH_MAX_DEPTH && n > BVH_MAX_NUM_OBJ)
	{
		float cost = alot;
		Objects objs1;
		Objects objs2;

		float tempcost;
		Objects temp1, temp2;

		// find the best split by checking along each axis
		for (int i=0; i<3; ++i)
		{
			// try 10 different splitting planes
			split_box(objs, i, 10, &tempcost, &temp1, &temp2);
			// if it's the cheapest so far, save it
			if (tempcost < cost)
			{
				cost = tempcost;
				objs1 = temp1;
				objs2 = temp2;
			}
		}

		// create child BVHs
		BVH *child1 = new BVH(m_depth+1);
		child1->build(&objs1);
		m_objects->push_back(child1);

		BVH *child2 = new BVH(m_depth+1);
		child2->build(&objs2);
		m_objects->push_back(child2);
	}
	else // this is a leaf box, so don't split any further
	    *m_objects = *objs;
}


// Find the best (lowest cost) splitting plane along a given dimension, using tightly 
// fitting candidate child boxes and surface area heuristics to compute cost.
bool BVH::split_box(Objects *objs, int axis, int nbins, float *cost,
				Objects *objs1, Objects *objs2)
{
	// Bin the objects along an axis (with one pass).
	// For each candidate split merge bins.
	// ex: a split between bins 3 and 4
	// child1 gets bin1+bin2+bin3, child2 gets bin4+bin5+...

	Vector3 range = m_max - m_min;
	float binsize = range[axis] / nbins;
	std::vector<Objects> bins;
	bins.resize(nbins);

	// divide objects into bins
	Vector3 bmin, bmax;
	int a=0, b=0;
	int n = objs->size();
	for (int i=0; i<n; ++i)
	{
		(*objs)[i]->getBounds(&bmin, &bmax);
		for (int j=0; j<nbins; ++j)
		{
			// TODO: we don't really need to do the last check, since whatever's left over should
			// be added to the last bin, but I'll deal with that later
			if ((bmin[axis] + bmax[axis])/2 < m_min[axis] + (j+1)*binsize + epsilon)
			{
				bins[j].push_back((*objs)[i]);
				break;
			}
		}
	}
	
	// examine potential BVHs

	objs1->clear();
	objs2->clear();
	*cost = alot;
	
	Objects temp1;
	Objects temp2;
	float tempcost;
	
	for (int i=0; i<nbins-1; ++i)
	{
		temp1.clear();
		temp2.clear();

		int j=0;
		for (; j<=i; ++j) // before the dividing plane goes in temp1
			temp1.insert(temp1.end(), bins[j].begin(), bins[j].end());
		for (; j<nbins; ++j) // and after goes in temp2
			temp2.insert(temp2.end(), bins[j].begin(), bins[j].end());

		// compute cost
		tempcost = BVH::computeCost(&temp1) + BVH::computeCost(&temp2);

		// if it's the cheapest configuration so far, save it
		if (tempcost < *cost)
		{
			*cost = tempcost;
			*objs1 = temp1;
			*objs2 = temp2;
		}
	}

	return true;
}

// compute the axis-aligned bounding box of a set of objects
void BVH::getBounds(Objects * objs, Vector3 *bvmin, Vector3 *bvmax)
{
	*bvmin = Vector3(alot);
	*bvmax = Vector3(-alot);

	Vector3 bmin, bmax;
	int n = objs->size();
	for (int i=0; i<n; ++i)
	{
		(*objs)[i]->getBounds(&bmin, &bmax);
		bvmin->x = std::min(bvmin->x, bmin.x);
		bvmin->y = std::min(bvmin->y, bmin.y);
		bvmin->z = std::min(bvmin->z, bmin.z);
		bvmax->x = std::max(bvmax->x, bmax.x);
		bvmax->y = std::max(bvmax->y, bmax.y);
		bvmax->z = std::max(bvmax->z, bmax.z);
	}
}


// Compute the cost of intersecting a bounding box around a given set of objects,
// using surface area heuristics.
// TODO: also factor in cost of boxes and cost of triangles (write Object.cost)
float BVH::computeCost(Objects * objs)
{
	// compute surface area of the bounding volume
	Vector3 bvmin, bvmax;
	getBounds(objs, &bvmin, &bvmax);
	Vector3 range = bvmax - bvmin;
	float sa = 2*(range.x*range.y + range.y*range.z + range.z*range.x);

	int n = objs->size();

	return sa * n;
}

// Intersect the BVH's outermost bounding box
// Assumes m_min and m_max in the bbox are correctly set to the min and max of the bounding box
bool BVH::bbox_intersect(const Ray& ray) const 
{
#ifdef RAY_STATS
	// accumulate global total
	g_stats.m_totalBBoxIntersections++;
#endif

#ifdef SSE
	__m128 PosInf = _mm_loadu_ps(INF4);
	__m128 NegInf = _mm_loadu_ps(N_INF4);
	__m128 o = _mm_loadu_ps( &ray.o.x );
	__m128 d = _mm_loadu_ps( &ray.d.x );
	__m128 boxmin = _mm_loadu_ps( &m_min.x );
	__m128 boxmax = _mm_loadu_ps( &m_max.x );
	
	__m128 t1 = _mm_mul_ps(_mm_sub_ps(boxmin, o), ray.invd);
	__m128 t2 = _mm_mul_ps(_mm_sub_ps(boxmax, o), ray.invd);

	__m128 tmin = _mm_min_ps(_mm_max_ps(t1,NegInf),_mm_max_ps(t2,NegInf));
	__m128 tmax = _mm_max_ps(_mm_min_ps(t1,PosInf), _mm_min_ps(t2,PosInf));

	__m128 tminy = _mm_shuffle_ps(tmin, tmin, 0x39);
	__m128 tmaxy = _mm_shuffle_ps(tmax, tmax, 0x39);

	tmin = _mm_max_ss(tmin, tminy);
	tmax = _mm_min_ss(tmax, tmaxy);

	__m128 tminz = _mm_movehl_ps(tmin,tmin);
	__m128 tmaxz = _mm_movehl_ps(tmax,tmax);

	tmin = _mm_max_ss(tmin, tminz);
	tmax = _mm_min_ss(tmax, tmaxz);

	int tempp = _mm_comige_ss(tmax,tmin);
	const bool ret = (_mm_comige_ss(tmax, _mm_setzero_ps()) &_mm_comige_ss(tmax,tmin));

	return  ret;

#else
	float tx1 = (m_min.x - ray.o.x)*ray.inv_d.x;
	float ty1 = (m_min.y - ray.o.y)*ray.inv_d.y; 
	float tz1 = (m_min.z - ray.o.z)*ray.inv_d.z;

	float tx2 = (m_max.x - ray.o.x)*ray.inv_d.x;
	float ty2 = (m_max.y - ray.o.y)*ray.inv_d.y; 
	float tz2 = (m_max.z - ray.o.z)*ray.inv_d.z;

	if(tx1 > tx2)
		std::swap(tx1,tx2);
	if(ty1 > ty2)
		std::swap(ty1,ty2);
	if(tz1 > tz2)
		std::swap(tz1,tz2);
	
	bool ret1 = (tx1 < ty2 && tx1 < tz2) && (ty1 < tx2 && ty1 < tz2) && (tz1 < tx2 && tz1 < ty2);
	return (tx1 < ty2 && tx1 < tz2) && (ty1 < tx2 && ty1 < tz2) && (tz1 < tx2 && tz1 < ty2);
#endif
}


// Traverse the BVH (this is recursive)
bool BVH::intersect(HitInfo& minHit, const Ray& ray, float tMin, float tMax)
{
	if (bbox_intersect(ray))
	{
		bool hit = false;
		HitInfo tempMinHit;
		minHit.t = MIRO_TMAX;
    
		for (size_t i = 0; i < m_objects->size(); ++i)
		{
			if ((*m_objects)[i]->intersect(tempMinHit, ray, tMin, tMax))
			{
				hit = true;
				if (tempMinHit.t < minHit.t)
					minHit = tempMinHit;
			}
		}
		return hit;
	}
	return false;
}
